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Abstract 

Numerical simulations of thin sheets undergoing large deformations are computationally challeng¬ 
ing. Depending on the scenario, they may spontaneously buckle, wrinkle, fold, or crumple. Nature’s 
thin tissues often experience significant anisotropic growth, which can act as the driving force for 
such instabilities. We use a recently developed finite element model to simulate the rich variety of 
nonlinear responses of Kirchhoff-Love sheets. The model uses subdivision surface shape functions in 
order to guarantee convergence of the method, and to allow a finite element description of anisotrop- 
ically growing sheets in the classical Rayleigh-Ritz formalism. We illustrate the great potential in 
this approach by simulating the inflation of airbags, the buckling of a stretched cylinder, as well as 
the formation and scaling of wrinkles at free boundaries of growing sheets. Finally, we compare the 
folding of spatially confined sheets subject to growth and shrinking confinement to find that the two 
processes are equivalent. 


1 Introduction 

Thin sheets are omnipresent in nature, technology and everyday life, appearing at virtually all length 
scales. Being much thinner in one than in the other two dimensions, they can develop an unparalleled, rich 
variety of deformation modes when subjected to external forces, spatial constraints, or intrinsic growth. 
They buckle, wrinkle, fold, and crumple. Numerical simulations are an often-indispensable approach for 
studying the complex interplay of these modes. The folding and crumpling of a piece of paper [1-7] 
and metal sheets wrinkling and crumpling in vehicle collisions [8-11] are two out of many examples. 
For many such problems, the finite element method (FEM) has shown to be amongst the most efficient 
and flexible tools, especially in cases with strong material nonlinearity, complex geometry, or anisotropy. 
Even though the Kirchhoff-Love theory [12] provides a simple kinematic description, numerically sound 
finite element implementations of thin sheets have turned out to be difficult and cumbersome in the past. 
These problems can be successfully overcome since the subdivision surface paradigm was introduced to 
the FEM [13,14]. 

Large deformations of soft thin tissue such as insect wings, plant leaves, cell membranes, or flowers are 
often induced by growth (or shrinkage) [15,16], inevitably leading to the development of residual stress 
[17,18]. In this paper, we present an extension of the Kirchhoff-Love theory to allow for anisotropic in¬ 
plane growth, which we implement with Loop subdivision shape functions. The combination of these two 
concepts grants access to a very straightforward and highly efficient, yet powerful and flexible numerical 
tool for the simulation of nonlinear thin sheet mechanics. Our approach accounts for the change of 
reference curvature when the surface grows, generalizing recently developed tethered mass-spring models 
[19,20]. The next section summarizes the mentioned extension. It is followed in the subsequent sections 
by a series of thin sheet problems that we solve using the developed FEM implementation. Special 
attention is paid to the formation of self-similar wrinkles along a plastically stretched free edge as well 
as on the scaling of single-wavelength wrinkles of growing cylinders similar to flowers. 
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2 The Kirchhoff—Love Sheet with Anisotropic Growth 

Let C E 3 be the stress-free undeformed (“reference”) middle surface of a sheet with small thickness h. 
Under the action of external forces or growth, the sheet deforms into a new configuration with middle 
surface U C E 3 . In the following, let Greek indices a,/?, 7 , 5 G {1,2}, and Latin indices i, j G {1,2,3}. 
Lower (upper) indices will denote covariant (contravariant) components. Moreover, let { 0 1 ,# 2 ,# 3 } be a 
curvilinear coordinate system, and let x( 6 d, 0 2 ) and x(# x , 0 2 ) be parametrizations of U and U, respectively 
(see Fig. 1). The material points p and p = x(p) in the reference and deformed sheet are parametrized 
as 

p(0\ 6> 2 , 0 3 ) = x(<9\ 6 2 ) + 6> 3 a 3 (<9\ 6 2 ) and p (<9\ 6> 2 , 6 3 ) = x(0\ 6 2 ) + 6> 3 a 3 (<9 1 ,6> 2 ), (1) 

where 0 3 G [— h/2,h/2]. y is a diffeomorphism that maps from the reference to the deformed material 
positions. The tangent spaces of U and U are spanned by the respective vector fields 

Fjnr f)~yc 

a a (9 1 ,9 2 ) sx„ = — and a a (9\9 2 ) = x, Q = —. ( 2 ) 

By virtue of the Kirchhoff assumption, straight material lines normal to the middle surface retain these 
properties as well as their length. They are determined by the unit normal vectors 


ai x a 2 . ai x a 2 

a 3 = —- —r and a 3 = - - 

ai x a 2 ai x a 2 


(3) 


X 



Figure 1: Reference, grown and deformed configurations of the sheet’s middle surface. 

The covariant components of the first fundamental forms follow as 

Uaf3 = a o: ■ and cta/3 = ■ a^, (4) 

while those of the second fundamental forms are given by 

baf3 — a 3 • a a ^/3 and — a 3 • a^^. (5) 

Assuming that the thin sheet obeys the St. Venant-Kirchhoff law of linear elasticity, the connection 
between its kinematics and energetics is provided by the Koiter energy density functional [ 21 , 22 ]. Let 
the sheet be characterized by Young’s modulus E and Poisson’s ratio v. The elastic energy U e of the 
Koiter sheet is obtained by integrating the energy per unit area over the middle surface: 

C4[x, X] = i J_ EE (H a ^ 5 a a0 a l5 + ^H a ^ s /3 a0 ^ dU, (6) 

where dU = |ai x a 2 | d 6 > 1 d 6 >2 . The Einstein summation applies to repeated indices. H is often referred 
to as the “elastic tensor”, and is given component-wise by 

H a ^ s = va aP a? s + E^(a a ^ 5 + a aS a^). (7) 

a = (a —a)/ 2 and /3 = b — b are the in-plane ( 2 x 2 ) membrane and bending strain tensors, respectively. 
The Koiter shell ( 6 ) can be extended to incorporate anisotropic growth through the multiplicative de¬ 
composition of the geometric deformation gradient Vy = F e F g [23,24] into a growth tensor F g and a 
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purely elastic response F e , that ensures continuity and compatibility of the body. Owing to the Kirchhoff 
constraints, we may write 


G 0 
0 T 1 


(G G R 2x2 symmetric), 


and the growth-modified elastic strains then simply read [25] 


a = - (G T aG 1 — a) , 
f3 = b — G _T 6G _1 . 


( 8 ) 

(9) 

( 10 ) 


We further augment the elastic energy (6) with an inertial term to capture the dynamics of the thin 
sheet. The kinetic energy reads 

[/k[x,x] = i [_hpx ■ xdfi, (11) 

1 Jn 

where p is the mass density of the sheet, and x = dx/dt is the velocity. Our aim is to find the minimizer 
x of the total energy U = U e + for given growth tensors F g or external driving forces. 


3 Finite Element Implementation with Subdivision Surfaces 

To account for out-of-plane bending rigidity, the bending term in Eq. (6) integrates the Gaussian and 
mean curvatures, which comprise second derivatives of the displacement field u = x — x, over the middle 
surface. For boundedness of the integral in the weak formulation, continuously differentiable finite element 
shape functions (G 1 -continuity) are needed. This requirement has proven very challenging in the history 
of shell finite elements, until Cirak et al. [13,14] have introduced Loop subdivision surfaces to the FEM. 
A fundamental difference to traditional finite elements is that subdivision surfaces gain C 1 -continuity at 
the expense of a larger support of the individual shape functions. Details on their implementation are 
given in Refs. [13,25]. 

Aside from guaranteeing convergence, subdivision surfaces allow a classical Rayleigh-Ritz formulation 
of the sought finite element deformation: no rotational variables are needed and the only unknowns 
are the three nodal displacements. Moreover, a single quadrature point per element is sufficient [13, 
14,25], rendering this finite element approach computationally highly efficient and flexible. Of course, 
increasing the number of quadrature points may assist in resolving strongly anisotropic growth fields or 
constitutive relations. We employ Loop subdivision surface shape functions here to minimize the total 
energy V numerically by solving Newton’s equation of motion with a standard predictor-corrector scheme. 
Subcritical viscous damping is added for numerical stability and equilibration. 


4 Inflated Pillows and Airbags 

As a first instance of folding and wrinkling, we reproduce the inflation of pillows and airbags from 
Ref. [14]. A square sheet with diagonal length d = 120 cm and thickness h = 1 mm and a circular 
sheet with radius R = 35 cm and thickness h = 0.4 mm are instantaneously pressurized with 5 kPa to 
buckle out of their flat initial configuration. The elastic moduli are given by E = 588 MPa, v = 0.4 and 
E = 60 MPa, v = 0.3, respectively. We exploit the reflection symmetry by simulating only the upper 
half of the geometry. The quasi-static equilibrium configurations are shown in Fig. 2, and a movie of the 
dynamic deformation is provided in the supplementary material. The square pillow features a distinct 
folding pattern in the middle of the four edges, resulting from the non-uniform distribution of Gaussian 
curvature: The edges are closer to the point of maximal uplift than the corners, thus getting curved more 
and pulled towards the center. They are laterally stretched and longitudinally compressed, yielding the 
observed folds. The circular airbag, on the other hand, is initially axisymmetric and therefore behaves 
differently. Wrinkles develop similarly to elastic plates stamped into curved cavities [26] or ultrathin films 
placed on fluid drops [27]. On top of these low-amplitude wrinkles, axisymmetry of the stress field is 
broken and large crumples occur that localize the geometrically imposed Gaussian curvature. 
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Figure 2: Inflated square pillow (a) and circular airbag (b). 


5 Buckling of a Stretched Cylinder 

A frequently studied buckling problem is the laterally stretched open cylinder [28,29]. Two opposite 
point forces of equal increasing magnitude F cause the cylindrical sheet with radius R = 4.953 cm, 
length L = 10.35 cm, thickness h = 0.94 mm, and free edges to first bend before snapping through at 
F c = 11.836 Eh 3 /R to the post-buckling regime, where further deformations are dominated by stretching. 
The elastic moduli are fixed to E = 10.5 MPa, v = 5/16. In Fig. 3, we plot the radial displacements of 
the points A, B and C, which are degenerate at the buckling threshold F c . The corresponding movie can 
be found in the supplementary material. 




Figure 3: (a) Reference state of the cylindrical sheet, (b) Equilibrium solution at F — 1.92F C . The 
stretching energy density is shown on a logarithmic color scale, (c) Normalized radial displacements 
vs. the rescaled point force. 


6 Boundary Instabilities and Wrinkling 

An interesting feature observed in plant growth is the occurrence of self-similar wrinkles along free tissue 
boundaries, such as the edges of flowers and leaves [19,20,30-33]. The morphology is very similar to the 
shape of torn plastic sheets, and apparently, both phenomena are characterized by a plastic longitudinal 
metric profile gi(z) = 1/(1 + z/l\ where / > 0 is a characteristic length scale and z > 0 is the coordinate 
perpendicular to the growing edge. Audoly and Boudaud [34] were able to show that the solution of the 
Foppl-von Karman equations on the edge of a free rectangular sheet with such growth profiles consists 
of self-similar wrinkles governed by odd integral scaling factors. The Foppl-von Karman equations 
are, however, geometrically limited as they don’t allow reentrancy. The present growing Koiter shell 
model allows us to numerically solve the problem with its full geometric nonlinearity taken into account. 
Consider a flat rectangular sheet of thickness h = 10 -4 , length L = 4, and width W = 1, initially lying in 
the xz plane. We clamp the long edge at z = W and constrain the short edges to stay at x = —L/2 , L/2 , 
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leaving them free to move in other directions. Plastic growth is imposed by setting the growth tensor to 

F g = diag(l + gi(z), 1, l) (12) 

in Cartesian coordinates (x,y,z), and we choose a characteristic length l = 40 h for the growth field in 
this example. Fig. 4 shows the resulting equilibrated configuration after growth (or tearing), and a movie 
showing the equilibration is provided in the supplementary material. The self-similarity of the free edge 
at £ = 0 is apparent, clearly resembling the wrinkling cascades observed in experiments [20,30]. 






Figure 4: Self-similar sheet boundary after growing according to Eq. (12). (a) Projection of the grown 
edge onto the xy plane, (b) 60° Fibonacci word fractal, (c) Koch snowflake. 


We have measured the fractal dimension of the grown edge depicted in Fig. 4(a) using the box counting 
method [35] and the self-similarity method [36]. In the former, the curve length L- in contained in a cubic 
box is determined as a function of the edge length Lbox of the box. A fractal curve is expected to 
scale as L[ n ~ L^ o f x . Such scaling is indeed observed with fractal dimension D$ = 1.15(1) (Fig. 5(a)). 
The scaling breaks down due to the influence of the clamped opposite edge, which introduces a global 
straightening effect when the box size is large (Lbox ~ W ). The second method is more robust to global 
orientation and is thus better suited here. The length L s of a piecewise linear path along the curve with 
segment size s is measured and expected to scale as L s ~ «s 1_jD . We find a self-similarity dimension 
D = 1.196(5) (Fig. 5(b)), which is very close to the Hausdorff dimension of a 60-degree Fibonacci word 
fractal (Fig. 4(b)), Dr = 1.2083 [37], and a bit lower than that of a triadic Koch curve (Fig. 4(c)), 
D h = 1.2619 [38]. 


^ (a) Box counting method 


o 



(b) Self-similarity method 



Figure 5: Fractal dimension of the edge at 2 = 0 of a thin sheet grown according to Eq. (12), measured 
with two standard methods. 


The metric profile used above is not the only one yielding wrinkled edges. When it comes to wavy 
flowers like certain orchids for instance, single-wavelength undulations instead of self-similar edges are 
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not uncommon. The feature causing wrinkle cascades is the presence of a non-constant geometric length 
scale defined by l geo {z ) = —g(z)/g'(z) [33]. The following families of growth fields will thus produce 
similar boundary instabilities: 


9l, P ( z ) x ( 1 + ^l) > 

oc (l-i) , 

On the other hand, an exponential growth field 

gi(z) oc exp (-|) , 


> 0, p > 0, z > 0), 

(13) 

> 0, p > 1, 0 < z <pl). 

(14) 

V 

0 

0 

IV 

2 

(15) 


yields only a single wavelength [19] because l geo (z) = l in this case. On some flowers, these undulations 
may be forced to integral wavenumbers n by angular periodicity of a single petal. Let’s hence significantly 
increase the characteristic length l and thickness h such that only a single wavelength A prevails even for 
growth in the form of Eqs. (13,14), and let’s consider a cylinder with height H and radius R instead of a 
flat plate. This change in geometry delivers dramatic consequences: A thin cylindrical sheet growing in 
circumferential direction according to g(z), where z is the cylinder axis, only breaks its axisymmetry if 
growth leads to a circumference that changes faster than the sheet’s metric can account for [31]. (Note 
that the excluded linear case p = 1 of Eq. (14) produces the excess cone [16, 39] in the limit R -> 0, 
which doesn’t wrinkle at the boundary because g[\(z) = 0.) A direct consequence of the Gauss-Bonnet 
theorem is that the axisymmetry is preserved as long as 




< 1, 0 <z<H, 


(16) 


and broken otherwise. The origin for this instability is the (non)existence of embeddings of the surface: 
According to Gauss’s Theorema Egregium, the creation or elimination of Gaussian curvature must be 
accompanied by in-plane stretching, which is traded for out-of-plane buckling if the sheet is sufficiently 
thin (see Fig. 6). 



Figure 6: Boundary instability of a circumferentially growing cylinder, (a) As long as Ineq. (16) holds, 
axisymmetry is preserved, (b) Same growth field g(z) oc exp(— z/l) as in (a), but larger prefactor. The 
axisymmetry is spontaneously changed to n-fold rotational symmetry C n ( h/R = 0.02, l/R = 1, n = 8). 
The rescaled bending energy density is shown in color, revealing the line where Ineq. (16) holds equally, 
(c) A relatively short cylinder (small H/l) with free boundaries also buckles away from the wrinkled edge, 
breaking C n symmetry further to C 2 . 


How does the number n of boundary waves scale with l,h,R when axisymmetry is broken, and is it 
universal for all positive, monotonically decreasing and strictly convex growth profiles satisfying 


lim /geo(^) — ^ ? 


z^O 


(17) 


Since the preferred wavelength A is a local feature independent of global geometry and topology (inde¬ 
pendent of i?), we may use the ansatz [33] 


A ^ 


h a l 


1—a 
geo ) 



(18) 
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On the other hand, geometry implies that 


A ~ —(1 + 0(0)). (19) 

since z = 0 is where Ineq. (16) is violated first, given that g' < 0 and g" > 0. After combining Eqs. (16-19), 
one thus finds a scaling for the number of wrinkles 


n ~ 




( 20 ) 


Up to the first term, which accounts for the mean curvature of the cylinder, this coincides with the scaling 
law reported in [34] for the wrinkling hierarchies in initially flat sheets, where a = 2/5 is found for the 
family gi(z) oc 1/(1 -\-z/l). Our numerical data, which best fits Eq. (20) with a = 0.39(2), shows that the 
wrinkles studied here fall in the same category, see Fig. 7. Moreover, the data collapse of all employed 
growth profiles on a single line indicates that the scaling is universal in this respect. For this numerical 
study, we used 

F g = diag(l, 1 + g(z), l) (21) 

in cylindrical coordinates (r,cp,z) with various growth profiles proportional to Eqs. (13-15), and with 
slowly increasing proportionality prefactors. A movie showing such spontaneous wrinkling for different n 
is provided in the supplementary material. 



Figure 7: Scaling of the wrinkle number n along the edge of a circumferentially growing cylinder. The 
mode with the lowest energy is the integer n closest to the power law (20), but excited modes also 
randomly occur and are metastable (data points lying significantly above or below the straight line). 


7 Confined Growth and Crumpling 

Many numerical simulations of thin sheets getting folded and crumpled inside of shrinking hollow spheres 
have been carried out recently [4-6,40]. An important result is that the high bulk stiffness of crumpled 
sheets is due to a network of vertices and lines carrying large mean curvature and bending energy [1]. 
But what if instead of externally forced compression, the thin sheet intrinsically grows inside of a fixed 
spatial container? Are the processes that crumple a plant leaf or petal growing inside a bud the same 
as for a piece of foil crumpled by hand? Indeed they are in the elastic limit, as the following simulation 
demonstrates. 

A thin circular sheet (radius R) is placed inside of a spherical cavity (radius R ). In the first setup, 
the container is shrunk, folding and crumpling the sheet into a ball of the size of the container. In the 
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second setup, the container sustains its size while the sheet undergoes uniform isotropic growth, both 
in plane and in thickness. The only important parameter for this problem is the Foppl-von Karman 
number 7 ^ (R/h ) 2 = 10 4 . Equivalent time scales are obtained by shrinking the sphere according to 
R(t) = R/ ( 1 +#(£)), where g(t) = At is the growth factor of the growing sheet. The growth rate A is chosen 
small enough to keep inertial effects negligible. We add repulsive contact forces penalizing volumetric 
overlap between any two pieces of the sheet. Initially, both sheets buckle to form a developable cone 
with a single vertex (see Fig. 8 , first column) that starts to nucleate at R/R ~ 0.53. The emerging 
ridge network of focused mean curvature (third column) and bending energy (fourth column) is the 
same in both scenarios. Compared to similar measurements [5], the cross correlation r = 0.89 of the 
mean curvature ridge patterns is very high. The differences are only local and of the order of the mesh 
resolution. A movie showing the crumpling process is contained in the supplementary material. 


R/R = 5/3 


R/R = 9/4 



Figure 8 : Comparison between shrinking confinement (bottom row) and growth in static confinement 
(top row), showing that the two processes are equivalent. k\ and k 2 are the principal curvatures of the 
middle surface. 


8 Conclusions 

Simulating the plastic growth of thin sheets can be numerically demanding. The finite element method, 
being perhaps the best-suited technique for anisotropies, lacked an expedient C 1 -continuous discretiza¬ 
tion, which is indispensable from a theoretical viewpoint, until subdivision surfaces were ported to it. We 
have extended the Kirchhoff-Love theory by arbitrary volumetric in-plane growth [25] and used the Loop 
subdivision surface paradigm to build a highly flexible and efficient numerical tool for the simulation 
of nonlinear thin sheet mechanics. Requiring no rotational variables, a thin sheet representation of this 
kind is remarkably simple to implement and superior to traditional approaches in terms of computational 
costs. A series of example simulations were carried out to demonstrate these strengths. In particular, 
we have quantified the self-similarity and scaling of wrinkles along plastically stretched free edges of thin 
sheets, and found that the problem of a growing confined sheet is equivalent to the narrowing confinement 
of a static sheet. 
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